Modified Kubelka-Munk equations for localized waves inside a layered medium 
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We present a pair of coupled partial differential equations to describe the evolution of the average 
total intensity and intensity flux of a wavefield inside a randomly layered medium. These equations 
represent a modification of the Kubelka-Munk equations, or radiative transfer. Our modification 
accounts for wave interference (e.g., localization), which is neglected in radiative transfer. We 
numerically solve the modified Kubelka-Munk equations and compare the results to radiative transfer 
as well as to simulations of the wave equation with randomly located thin layers. 
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I. INTRODUCTION 

The most basic mesoscopic theory that attempts to 
explain multiply-scattered wave energy is the theory of 
radiative transfer (RT). In a ID layered medium, RT is 
equivalent to the well-known Kubelka-Munk (KM) equa- 
tions P, Q because the KM equations are in essence a 
two-flux theory and in ID there are only two directions 
(up and down). However, due to the inevitability of wave 
interference in ID Q , RT is unable to accurately predict 
all aspects of energy transport in randomly layered media 
[H . Wave interference is explicitly ignored in RT 0, [f| 
and leads to the phenomenon of wave localization, as de- 
scribed by many authors previously 0, H, 0, HI- Wave 
localization is of primary importance for topics such as 
the interaction of electrons with disorder 0] (e.g, the 
metal-to- insulator transition), the transmission of light 
through randomly layered structures (such as a stack of 
transparencies |7[), and the late-time behavior of seismic 
recordings at volcanoes [Toj . 

Previous studies have been devoted to understanding 
RT in layered media, despite its neglect of wave interfer- 
ence. Hemmer [ll[ may have been the first to solve for the 
Green's function of RT in ID, as pointed out by Paass- 
chens [l2j]. The application of ID RT to vertical seismic 
profiles has been the subject of work by Wu [H, Q and 
Wu and Xie fl5ll. Sato and Fehler [l6| have discussed ID 
RT and Sato [17| has derived the solution of the Green's 
function of RT in ID using the integral form, instead of 
the differential form used by Hemmer [Til ] . Building upon 



the work of Wu and Xie [15[ , which centered on station- 



ary RT in layered media, the time dependent case has 
been recently considered [HI . Bakut et al. have gen- 
eralized the Green's function of ID RT in homogeneous 
media to the case of a medium composed of piecewise ho- 
mogeneous layers. Though ID RT has been thoroughly 
understood in the course of these studies, how wave in- 
terference changes the picture - from the point of view 
of RT - has so far not been covered. It is the aim of the 
present work to properly incorporate wave interference, 
and hence the phenomenon of wave localization, within 



the framework of RT for the case of layered media. 

Wave localization in ID systems has received consider- 
able attention, both theoretically and experimentally. As 
a result, several different techniques have been applied. 
Among the most widely used is random matrix theory 
coupled with Furstenberg's theorem @, [M HH, HH . This 
approach deals with the wavefield itself for a single real- 
ization of randomness by using so-called "self-averaging" 
quantities [23|. Given an ensemble of random realiza- 
tions, these quantities converge (closely) to their aver- 
age in a single realization, provided the realization in- 
cludes enough scatterers. In spite of its ability to model 
the wavefield itself, random matrix theory is basically a 
stationary theory and it relies on a limiting procedure, 
Furstenberg's theorem, which takes the limit as the num- 
ber of matrix products (i.e., scatterers) becomes infinite. 
Furthermore, random matrix theory is primarily limited 
to ID systems. Historically, this fact has led to a discon- 
nect in the prevailing theoretical treatment of multiple 
scattering in ID (random matrix theory) versus 2D and 
3D (RT). 

Significant progress has been made recently toward in- 
corporating wave interference into RT (at least within the 
diffusion approximation) using the self-consistent (SC) 
theory of Anderson localization. In fact, a ID version 
of the SC-theory has been studied analytically [2J] . The 
SC theory is different from random matrix theory in that 
it predicts the late time spatial and temporal evolution 
of the mean wavefield intensity (the squared wavefield) 
for an ensemble of random realizations. The crux of the 
SC theory is that it attempts to include the effects of 
wave interference by using a so-called "self-consistent" 
expression for the diffusion constant, an idea originally 
popularized by Vollhardt and Wolfle [9(. 

Here, we attempt to include the effects of wave in- 
terference by deriving the ID RT equations from a fun- 
damental level, using a procedure first demonstrated by 
Goedecke @j. We find that once properly modified, the 
ID RT equations (also known as the KM equations Q) 
can account for interference effects such as wave local- 
ization. We call these new equations modified KM equa- 
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tions. Thus, we are able to correctly account for wave 
interference within the framework of RT, at least in ID. 
We also show that the predictions of the modified KM 
equations agree with predictions of random matrix the- 
ory, namely the expected exponential decay of the steady- 
state transmission coefficient with sample size. We fin- 
ish by testing and verifying the modified KM equations 
through a comparison with numerical simulations of the 
wave equation. In contrast to the ID version of the SC 
theory [2J] , the modified KM equations hold for all times 
and model both the total intensity and the intensity flux. 
At the end, we comment on the prospects of generalizing 
the ID theory to higher dimensions, especially 3D where 
the notorious and interesting transition from extended to 
localized wave propagation occurs. 



II. THE SCATTERING MATRIX WITH 
INTERFERENCE TERMS 

We aim to derive equations similar to the ID RT equa- 
tions, or KM equations, but with the explicit inclusion 
of wave interference. Although it is not necessary we 
assume in the following that there is no absorption for 
simplicity. For a layered medium made up of thin lay- 
ers, or ID scatterers, embedded in a homogeneous back- 
ground medium, the scattering matrix relating incident 
and scattered waves at scatterer n is 
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where r and t are the reflection and transmission coeffi- 
cients of a scatterer, S n +i and S n +i are the downward 
and upward propagating complex wave amplitudes at the 
base of scatterer n, and S_ n and S_ n are the downward and 
upward propagating complex wave amplitudes at the top 
of scatterer n, as shown in Figure [TJ Note that the com- 
plex wave amplitudes at the top of scatterer n, S_ n and 
S_ n , are related to the complex wave amplitudes at the 
base of scatterer n — 1, S n and S n , by simple phase ad- 
vance or delay 
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where V Z is the delay operator associated with the 
propagation time between scatterers n and n — 1 po] ]. 
From equation it follows that \S_ n \ 2 = \S n \ 2 and 
\S_ n \ 2 = |SVi| 2 . By taking the squared magnitude of the 
two equations making up the scattering matrix, equa- 
tion {TJ), and adding and subtracting them, we thus arrive 
at the equations 
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FIG. 1: The up- and down-going waves near scatterer n. We 
consider a random medium consisting of thin layers, or ID 
scatterers, of thickness d and local wavenumber k embedded 
in a homogeneous background medium of wavenumber fco. 
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For a ID scatterer embedded in a homogeneous medium, 
two identities exist: |r| 2 + \t\ 2 = 1, or conservation of 
energy, and r*t + t*r ~ 0, as shown by Ursin [2^|. With 
these identities, expressions ((3]) and (j4]) can be simplified 

as 
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Note that the interference terms are not present in equa- 
tion iJSJ. This equation states that the energy flux be- 
tween scatterers n and n — 1 equals that between scatter- 
ers n and The principle of energy flux conservation 
in layered media has been used by Claerbout to derive 
the method of "acoustic daylight imaging" [2^, in chapter 
8]. In contrast, equation ©, which describes the local 
change in total intensity on either side of a scatterer, 
contains an interference term. This term depends on the 
correlation between the downward propagating wavefield 
at the top of scatterer n {S_ n ) and the upward propagat- 
ing wavefield at the base of scatterer n (S n +i)- 

We continue by averaging equations (0) and Q over 
ensembles of randomly placed scatterers. We denote en- 
semble averages with brackets, () and thus the ensem- 
ble average of the squared magnitude of the down-going 
wavefield between scatterers n— 1 and n by {\S n \ ) = In 
and so on for other wavefield quantities. With ensemble 
averaging, we obtain from equations ([5]) and ([6]) that 
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The stationary ID RT equations result from these equa- 
tions by first assuming zero correlation in phase between 
wavefields propagating in opposite directions at scatterer 
n i {B. n ^n+i) = 0, followed by taking a limiting procedure 
to move from discrete to continuous variables, as shown 
by Goedecke 0. 

It is well known, however, that wave interference causes 
the term (S_ n S* l+1 ) to be nonzero, especially in ID. We 
thus include the term containing (S_ n S* +1 ), and therefore 
account for interference effects, by considering the direc- 
tional wavefields on either side of a planar source within 
a ID random medium. We take the depth axis (z-axis) 
positive downward. First, consider the situation above 
the source depth z s . There, the up-going wavefield S_ n 
is the wavefield incident from the source direction and it 
can be related to the down-going wavefield S_ n using the 
reflection coefficient i? t for the entire random medium 
above scatterer n as 12711: 



(9) 



Moreover, the up-going wavefield S n +i can be related to 
the up-going wavefield on the other side of scatterer n, 
S_ n , using the reflection and transmission coefficients of 
a single scatterer, r and t, and the reflection coefficient 
for the entire random medium below scatterer n, denoted 
R 2 , as: 
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Equations ([9j) and JTDJ) relate the wavefields S n +\ and 
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Substituting this relationship for S_ n into equation ([5]) 
and assuming the ensemble averaging can be distributed 
as follows: 

Im(-^-|S n+1 | 2 ) = (|5„ +1 | 2 )Im(-^-), (12) 
1 — rR% 1 — rR.2 

equation j8]) above the source becomes 
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Applying the same considerations to the situation be- 
low the source depth z s means that the direction of the 
wavefield incident from the source is the opposite of the 



case just shown. In addition, the roles of the terms Ri 
and i?2 are different: Ri is now the reflection coefficient 
for the entire random medium beneath scatterer n and i?2 
is the reflection coefficient for the entire random medium 
above scatterer n. This convention maintains the same 
relation between i?i and i?2 and the direction of the in- 
cident wavefield as was used previously. Thus, starting 
with RiS n +i — SVi+i, equation ((5|) below the source is 
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For a single realization of the ensemble, the R\ and i?2 
here are not necessarily equal to the R\ and R2 consid- 
ered previously. However, the ensemble averages of the 
reflection coefficients on either side of the source are the 
same since the spacings of the scatterers above and below 
the source are drawn from the same random distribution. 
From equations (|T3| and (fl4| . the two situations differ 
not only by the direction of the wavefield present in the 
last term (/„ or I n +i), but also by a sign change in the 
last term (since sgn[Im( )] = -sgn[Im( T ^-)]). 



III. MODIFIED KM THEORY: 
STATIONARY CASE 



THE 



With these two cases (above and below the source), 
we now take the limiting procedure - as discussed by 
Goedecke Q - to move from the discrete to the continu- 
ous case. We examine here the case below the source and 
then state the result for the case above the source, since 
the procedure for the two cases is the same. First, note 
that In+i and I n +i are defined at the base of scatterer 
n, just as /„ and /„ are defined at the base of scatterer 
n — 1. We define the average spacing between the scat- 
terers as p" 1 and thus the number of scatterers per unit 
depth is p (the number density). Multiplying both sides 
of equation (| 14[) by p results in 
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As pointed out by Goedecke Q, the term on the l.h.s. of 
equation (|15p becomes a spatial derivative when making 
the transition to a continuous depth variable np^ 1 —>■ 
z. Therefore, the directional wavefields become functions 
of z, that is, I n = Id{z) and I n = I u {z) where Id and 
I u are the ensemble-averaged down-going and up-going 
intensities. Therefore, equation (|15p becomes 
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We further simplify equation (fl6|) by denning the 
ensemble-averaged total intensity I t = Id + I u an d the 
ensemble-averaged intensity flux If = I d — I u , This sim- 
plifies equation (|TB|) as 



dlt „ M 2 T , T lm(r*t) T , mt* 



T - r*R* 



(17) 



We finally define the scattering mean free path £ s , the 
localization length £i oc , and a dimcnsionless parameter 
B as 



is P W 
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Using these parameters, we can rewrite equation (|17p 
concisely as 
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We have chosen the definitions in equation (fT5)) in a 
manner consistent with the usual definitions for these 
quantities 0, H|. F° r instance, regarding the quantity 
/ 9 l r | 2 /K| 2 j m the weak scattering limit (\t\ 2 sw 1) we find 
that 
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for z < z s 
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These equations comprise the modified KM equations 
in the stationary case. Equation (|24p may be rewritten 
more concisely as 
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where the quantity It — \If\ is either 2I U for z > z s or 2/^ 
for z < z s : it is twice the intensity propagating back to- 
ward the source. Equation (f2"S"|) shows that the inclusion 
of wave interference in the KM (or RT) equations leads 
to two additional terms which affect the average total in- 
tensity in different ways. The first term containing l/£i oc 
on the r.h.s. of equation ([23]) causes the coherent inten- 
sity to decay more rapidly than when wave interference 
is neglected. Furthermore, the second term containing 
1/^/oc on the r.h.s. causes the spatial distribution of the 
incoherent intensity to be entirely different than in the 
the case of no interference (as demonstrated later in a nu- 
merical example) . The form of equation (]25[) allows the 
identification of the extinction mean free path (the decay 
of the coherent intensity), l/l ex t — B/£ s + l/£i oc . This 
insight is possible since the quantity It — \ If | in equa- 
tion (|25)) is zero for the coherent intensity. Note that 
the coherent intensity decays exponentially even without 
interference effects (£i oc — > oo, or RT) due to scattering 
out of the forward direction. 



B = 



t-1 



(21) 



Thus, B is a dimcnsionless parameter describing the di- 
rectionality of the scattering For isotropic scatter- 
ing, B = 1/2. In addition we define 



1 



p(kl 2 + |t-i| 2 ), 



(22) 



consistent with what we know for the ID scattering cross 
section from Sheng [jj: <r s = |r| 2 + |i— 1| 2 . Therefore, 
from equation (|2"2"|) , we can identify the factor p{\r\ 2 + 
\t — 1| 2 ) = pa s . In the weak scattering limit, it is well 
known that £ s — l/pa s . Thus, our definition for £ s in 
equation ([22]) is consistent with the usual definition of £ s 
in the weak scattering limit. Appendix lAl demonstrates 
the consistency of the definition for £i oc as it appears in 
equation (fl9)l based on the relation in equation (fTSf . 

We have just shown how to apply the limiting proce- 
dure to equation fT4")) . Applying the same limiting pro- 
cedure to equations (J7J) and (fT3"|) gives all of the neces- 
sary stationary transport equations, which we summarize 
here: 



dlf_ 

dz 



0, 



(23) 



IV. MODIFIED KM THEORY: THE 
TIME-DEPENDENT CASE 

Having derived the modified KM equations for the sta- 
tionary case in equations ([23)) and (|2"3")l , we will turn 
our attention to the time-dependent (dynamic) case. 
Given the current coordinate system for z, this is accom- 
plished by noting that dl u /dz = dl u /dz + v~ 1 dl u /dt 
and dld/dz = dld/dz — v~ 1 dl d /dt, where v is the ve- 
locity of energy transport (the energy velocity) . Includ- 
ing the presence of planar isotropic (zero net down-going 
component) sources [llj], we obtain the following time- 
dependent equations 



dh 
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v dt 
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dif + ldit _ r 

dz v dt v' 
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^loc 
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(26) 



(It-\If\), 



(27) 

where T is the isotropic (omnidirectional) source term 
[HI]. Note that for £i oc — > oo (no wave interference), 
equations ([26]) and ([27]) are the same equations as have 
been studied previously by others within the context of 
RT in layered media 0, M, H [3 ■ 
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With equations (|26|) and (f27|) . which are the modified 
KM equations in the time-dependent case, we proceed 
to numerically solve the equations for two cases: with 
interference and without [li oc — * oo, or RT). These cases 
are compared to ensemble averages of simulations of the 
wave equation. 

V. NUMERICAL SIMULATIONS 

Our numerical solution of equations (|26p and (I27p ex- 
ploits staggered grid finite difference methods [28|]. In 
this technique, we calculate the total intensity I t and the 
intensity flux If on different spatial grids that have been 
shifted by half of a grid-spacing and on different tempo- 
ral grids shifted by half of a time step. Our purpose is to 
test whether the modified KM equations, equations ([26]) 
and (|27|) . predict the results of a wave-based simula- 
tion. Therefore, we also simulate the wave equation for 
normally-incident plane waves in a layered medium, ex- 
cited by a planar force source: 

1 d 2 u d 2 u 1 _ . . . .„„. 
^W-d^ = ^/ Wmz) > (28) 

where u is the displacement field as a function of time 
t and spatial coordinate z, c is the phase velocity, p is 
the density of the medium, F is a dimensionless constant 
related to the strength of the forcing function, w(t) is 
the dimensionless source time function, and S(z) is the 
spatial delta function. We simulate equation l|28p by the 
finite-difference method using centered, second-order ap- 
proximations for the derivatives. The details of the nu- 
merical implementation have been previously discussed 
in Haney et al. [l8j . 

The setup of our numerical simulation is as follows: 
for a single realization, we place 50 scatterers randomly 
over a depth range of L = 200 m. The scatterers are 
lower in propagation velocity (1000 m/s) than the back- 
ground medium (2000 m/s) but have the same density. 
We excite a source in the center of the 200-m range at 
z s = m. At the ends of the 200-m range are absorb- 
ing boundaries. To obtain ensemble averages of the total 
intensity, we first bandpass filter our numerical results 
from a single realization with a Gaussian filter peaked 
at 500 Hz. We filter in the frequency domain since the 
transport properties (i.e., £ s and £i oc ) are strongly de- 
pendent on frequency. In other words, equations (|26|) 
and (|27|) model the wave experiment in a particular fre- 
quency band. After filtering, we square the wavefield for 
each of 100 simulations with randomly placed scatterers 
and then add the squared wavefields. From the ensemble- 
averaged wavefields, we estimate the extinction mean free 
path £ ext from the decay of the coherent intensity and 
the localization length £i oc from the exponential decay of 
the incoherent intensity away from the source. We find 
i ext = 38.1 ± 0.5 m and £ loc = 57.2 ± 1.7 m. These two 
parameters enter into equations (|26p and (j2"7| . We fur- 
ther find that the energy velocity of the coherent wave is 



only slightly altered from the phase velocity of the back- 
ground medium (2000 m/s), which is expected since the 
scatterers we employ are ID versions of Rayleigh scat- 
terers (B = 1/2, with thickness d much less than the 
dominant wavelength) ;3|, and hence are not resonant 
scatterers. 

The thick black line in Figure is the total intensity 
from the numerical solution of the standard KM equa- 
tions (KT, £i oc — > oo), with the wave simulation shown 
as the thin blue line. Note that these snapshots are log- 
arithmic in intensity. Strong localization effects are ev- 
ident in the wave simulation as seen in the sharp expo- 
nential peak in the total intensity at the source position 
at later times. This behavior is not captured in the so- 
lution of the standard KM equations, which predict that 
the total intensity is flat around the source position. In 
addition, the standard KM equations significantly under- 
predict the decay of the coherent wave. The discrepancy 
between standard KM theory and the wave simulation is 
most evident at t = 0.11 s in Figure^d), where the wave 
simulation shows a concentration of total intensity near 
the source position. 

The simulation for the modified KM equations is the 
thick black line in Figure [3j In contrast to the standard 
KM equations (or RT), the modified KM equations cap- 
ture the exponentially-peaked behavior of the total inten- 
sity near the source position and, at all times, agree well 
with the wave simulation. Thus, the modified KM equa- 
tions are capable of modeling the transport of intensity 
in ID localized media, where interference effects cannot 
be ignored. It is worth emphasizing finally that both the 
standard KM solution in Figure [2] and the modified KM 
solution in Figure [3] satisfy global energy conservation. 



VI. CONCLUSION 

With a proper modification to the well known Kubelka- 
Munk (KM) equations, we are able to accurately describe 
the transport of wave intensity in a ID layered medium at 
all times, even when interference effects dominate (e.g., 
wave localization). This is confirmed by numerical sim- 
ulations comparing wave simulations and the modified 
Kubelka-Munk equations. In the future, we plan to ex- 
tend our approach, which currently uses only two fluxes, 
to a theory valid for 2D and 3D disordered media. One 
approach to this extension would utilize higher dimen- 
sional discrete flux theories as described by Cwilich [29| . 
Such a transport theory will be capable of simultaneously 
describing the propagating coherent intensity, the inten- 
sity flux, and the localization transition in 3D. 
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FIG. 2: Comparison of numerical results for ensemble-averaged wave propagation (thin blue line) and standard KM theory, or 
RT (thick black line). The various panels show: (a) t = s, (b) t — 0.02 s, (c) t = 0.045 s, (d) t — 0.11 s. The source time 
function is zero-phase and hence acausal (symmetric about t = s). 
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APPENDIX A: VERIFICATION OF THE 
LOCALIZATION LENGTH 

Here, we give credence to our use of the term localiza- 
tion length £i oc as it appears in equation p^|) . By doing 
so, we justify the expression for £/ oc in equation (fT8|) . We 
proceed by finding the stationary transmission coefficient 
for a slab geometry in the case of interference. We adopt 
the approach shown by van Rossum and Nieuwenhuizen 
[30| , wherein the authors derived the stationary transmis- 
sion coefficient T for the case of no interference. In that 
case, T(L) ss z e / L, where z e is the extrapolation length 
outside the slab and L is the thickness of the slab [30] . fn 



analogy to electronic systems, the behavior T(L) z e / L 
is an expression of Ohm's law. 

We begin by taking the stationary version of equa- 
tions (Hi) and (H3 



dh 

dz 



D 



ho 



dz 



r 

3 

V 



sgn(z - z s ) 

^■loc 



(It 



'/I 



(Al) 



(A2) 



where, as discussed before in reference to equations (j2"rj)l 
and (|2"T)) . r is the isotropic (omnidirectional) source term 
[l8j . Let a single stationary (planar) source act at depth 
z s , such that T = S(z — z s ). Based on physical con- 
siderations, we know that /( is symmetric (even) about 
z = z s and If is antisymmetric (odd). This, together 
with the fact that T = 8(z — z s ) in equation (|A1|) . leads 
to the relation // = sgn(z — z s )/2v and therefore that 
sgn(Jf) = sgn(z - z s ). Since sgn(//)|//| = If, equa- 
tion (|A2[) may be written as 



dh 

dz 



2B 

17 



i 



sgn(z - z s ) 

^loc 



(A3) 



Solving this equation for If allows a substitution for If 
in equation (|A1|) . At depths away from the source (for 
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z 7^ z s ), this gives an expression in terms of It only 

d It 1 d r , . i . , . 

^ri + ^izi s ^ z ~ z ^^°- ( A4 ) 

UfZ tloc "Z 

For the case of no interference, h oc — > oo, equation (|A4|) 
reduces to the Laplace equation (the diffusion equation 
in the stationary case). However, when interference is 
taken into account, the equation is a modified Laplace 
equation. 

In preparation for an application of the standard ap- 
proach shown by van Rossum and Nieuwenhuizen [30(, 
we proceed by investigating how the extrapolation length 
outside of a slab of thickness L changes when interfer- 
ence is accounted for. We use the well-known approach 
of Morse and Feshbach [3l| to define the extrapolation 
length. Consider a slab containing randomly located thin 
layers extending from z = to z = L. Outside of this in- 
terval, the medium is homogeneous. Take a planar source 
of intensity at some z s < 0, outside of the slab. Since the 
entire slab extends over z — to z = L, we have z s < z 
and thus sgn(z — z s ) = 1 for all points z within the slab. 
In this case, equations (|A3[) and (|A4[) are, for all points 
z inside of the slab, given by 



dh 

dz 



2B 

17 



i 



If 



IT" 

L loc 



and 



d 2 i t 

dz 2 



j_dh 

hoc dz 



(A5) 



(A6) 



At the far end of the slab z — L, we require there 
be no up-going intensity = (I t — If)/2 = 0. Now 
using equation (IA5j) . we substitute for If in the relation 
I t — If = 0, giving an equation in terms of It only 



It + 



hoc + It 



dh 

dz 



^loc 



0. 



(A7) 



where we use the transport mean free path £ tr = £ S /2B 
to make the notation concise [HI]. We can rewrite equa- 
tion E3 as 



l t + a^=0 

dz 



(A8) 



where a = £trhocl {hoc + 2£ tr )- Equation (|A8|) means 
that, near z = L, I t ~ C(L + a — z)/a where C is a 
dimensioning constant. Within this approximation, It = 
at z = L+a\ therefore, the extrapolation length z e - the 
distance outside of the slab where I t vanishes - is equal to 
a. That is, z e — itr^ioc/ {hoc + 2£ tr ). One can see in this 
expression that, for no interference {£i oc — > oo), z e — £ tr 
which is the usual extrapolation length encountered in 
ID when interference is neglected [l8| ]. 



At the side of the slab on which the source of intensity 
is incident, at z = 0, we require the down-going intensity 
to be equal to the incident intensity, Iq. Th us, Id = 
{It + If)/2 = Iq at z = 0. Using equation (|A5|) . we 
substitute for // in the relation {I t + If)/ 2 = I , giving 



£tr£, 



tr^loc 



dz £i oc 
Equation (|A9|) may be rewritten as 



= 2/n 



(A9) 



(A10) 



8I t _ 21 {£ loc + £ tr ) 
oz hoc 

Near z — 0, the solution is approximately given by 

2lo{£ l oc + £, r )\z + J (AU) 

Within this approximation, at a distance equal to the (in- 
terference adjusted) extrapolation length outside of the 
slab, z = — z e , It is therefore given by 

£loc{£loc + 2£ t r) ° ' ^ ' 

where we represent the term containing £i oc and £ tr by 
A. 

Following the method employed by van Rossum and 
Nieuwenhuizen [3fJ , we seek to solve equation (|A6|) with 
the boundary conditions I t — and I t — 2/ A at the 
(interference adjusted) extrapolation lengths, z = L + z e 
and z = —z e respectively. The solution is 



It{z) = 2I A 



,-z/tloc 



-(L+ Ze )/t lo 



e/tloa 



-(L+z e )/i lo 



(A13) 



The steady-state transmission coefficient, T{L) = I t {z 
L)/I Q , is 



T{L) = 2Ae- i/£ '" 



2Ae" 



e -(L+z c )/e loc 
1 - e - z '/ tl 



.(A14) 



where the approximation is for z e <C L. This expres- 
sion shows that when interference is taken into account 
the steady-state transmission coefficient goes down ex- 
ponentially as a function of the length of the slab L - a 
hallmark of localization in the stationary case. This be- 
havior is in stark contrast to T{L) w z e /L [30( obtained 
when interference is neglected {£i oc — ► oo)- The fact that 
the length scale controlling the exponential decay with L 
in equation (|A14[) is £i oc supports the use of this term in 
equation (|19[) and, as a consequence, the expression for 
£i oc in equation (TTSJ). 



